Modeling oscillatory Microtubule— Polymerization 
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Polymerization of microtubules is ubiquitous in biological cells and under certain conditions it 
becomes oscillatory in time. Here simple reaction models are analyzed that capture such oscilla- 
tions as well as the length distribution of microtubules. We assume reaction conditions that are 
stationary over many oscillation periods, and it is a Hopf bifurcation that leads to a persistent 
oscillatory microtubule polymerization in these models. Analytical expressions are derived for the 
threshold of the bifurcation and the oscillation frequency in terms of reaction rates as well as typical 
trends of their parameter dependence are presented. Both, a catastrophe rate that depends on the 
density of guanosine triphosphate (GTP) liganded tubulin dimers and a delay reaction, such as the 
depolymerization of shrinking microtubules or the decay of oligomers, support oscillations. For a 
tubulin dimer concentration below the threshold oscillatory microtubule polymerization occurs tran- 
siently on the route to a stationary state, as shown by numerical solutions of the model equations. 
Close to threshold a so-called amplitude equation is derived and it is shown that the bifurcation to 
microtubule oscillations is supercritical. 

PACS number(s): 47.54. -f-r, 87.10. -He, 87.16.Ka, 87.17.Aa 



I. INTRODUCTION 

Microtubules are cylindric filaments that are used in 
cells for many different purposes, being vitally involved 
in cell motility and division, in organelle transport, and 
in cell morphogenesis and organization [1]. The precise 
ways in which microtubules achieve their amazing vari- 
ety of cellular functions is not fully understood yet. Mi- 
crotubules in cells are generally dynamic, they assem- 
ble, disassemble or rearrange on a time scale of minutes. 
GTP (guanosine triphosphate) hydrolysis is apparently 
the driving force of microtubule physiology. 

The rich non-equilibrium dynamics of microtubules, 
including the nucleation and polymerization kinetics etc. 
[2,3] is attracting considerable attention, both experi- 
mentally and theoretically [4-18]. Two phenomena in 
this area, the dynamical instability of microtubules [4] 
and the oscillatory polymerization [5-12] challenge theo- 
retical modeling already for a while [13-18]. 

Oscillations during microtubule polymerization have 
been observed either when GTP is regenerated enzymat- 
ically from endogenous GDP (guanosine diphosphate) 
[5,7,9,11] or when some amount of GTP is provided at 
the beginning or during an experiment. In the latter case 
oscillations occur only as a transient, because GTP is ei- 
ther consumed or some reactions steps may be inhibited 
due to the accumulation of GDP [8] . If both possibilities 
are combined, the length of a transient regime depends 
on the initial concentrations of GTP and GDP and on 
the capacity to regenerate GTP. Present models for mi- 
crotubule polymerization focus mainly on a description 
of transiently occurring oscillations and the solutions of 
the respective models are mostly numerical [6,15,17]. 

In recent in vitro experiments, however, the capacity to 
regenerate GTP has been enhanced and extended up to 
several hours [19]. Compared to a typical oscillation pe- 



riod during microtubule polymerization, which is of the 
order of a minute, the reaction conditions in these ex- 
periments are almost quasi-stationary over a long range 
of time. Therefore we focus on modeling microtubule 
polymerization for time-independent regeneration condi- 
tions. As starting point we take common reduced models, 
where several elementary processes of the real biochemi- 
cal reaction are described by a few effective reaction steps 
as explained in Sec. II, cf. Refs. [6,7,15,17]. Reductions 
of complex chemical reaction schemes are quite common 
and a famous example is the so-called oregonator [20] 
that is a reduced model for the legendary Belousov- 
Zhabotinsky (BZ) reaction [21]. However, since micro- 
tubules are long filaments, there are essential differences 
between the polymerization of microtubules filaments 
and common chemical reactions. For instance micro- 
tubules may undergo an orientational ordering transition 
beyond a critical filament density [22], a phenomenon, 
which doesn't occur in common chemical reactions. Ac- 
cordingly, the length distribution of microtubules is ex- 
plicitly taken into account for all variants of models in- 
vestigated in this work. Such models are the basis of fu- 
ture work on interesting pattern-formation phenomena 
related to the interplay between orientational ordering of 
the filaments and the kinetics involved in the filament 
growth [23]. 

In addition, we focus on model variants that include 
the possibility of an oscillatory microtubule polymeriza- 
tion and that allow analytical approaches. However, the 
reaction steps, such as nucleation, growth and decay of 
microtubules or the rate limiting factors of oligomer de- 
cay or tubulin regeneration, which have been identified to 
be crucial for oscillations [5-8,11], are taken into account. 
Moreover we address the question whether microtubule 
oscillations occur transiently or in a persistent manner 
beyond a Hopf bifurcation. Whether such a Hopf bifur- 
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cation takes place super- or subcritically is investigated 
in terms of the so-called amplitude expansion. 

It is not a major goal of this work to achieve quanti- 
tative agreement between the results obtained with phe- 
nomenological models and experimental measurements. 
However, since the present understanding of the mecha- 
nism leading to oscillatory microtubule polymerization is 
incomplete, reduced models may be an appropriate tool 
for working out typical trends that may be testable in 
experiments. For comparison it is very helpful that for 
reduced models trends may be worked out analytically 
and may be presented by simple formulae. A number of 
spatiotemporal phenomena involving microtubule poly- 
merization call for a better understanding too [10,24,25], 
but also in this case simple and effective models are in- 
dispensable in order to keep the modeling tractable [26] . 

At the transition to oscillatory polymerization the sta- 
tionary state becomes sensitive against small perturba- 
tions, which grow or decay exponentially, oc e"^*. Here the 
exponential factor a = Ur^ioJc is the sum of the so called 
growth rate cTr and the oscillation frequency Wc 7^ 0. Be- 
low the bifurcation point the growth rate Ur- < is neg- 
ative and the perturbations arc damped. Beyond the bi- 
furcation point Ur is positive and the stationary polymer- 
ization state is unstable against oscillatory perturbations. 
Hence the Hopf bifurcation to oscillatory polymerization 
takes place when the real part ct^ of both roots passes 
zero. The investigation of the polymerization dynam- 
ics beyond the Hopf bifurcation requires in most cases a 
numerical analysis of the basic reaction equations. How- 
ever, close to threshold ar is small and the oscillation of 
the polymerization, described by the real part of 6*'^'=*, is 
much faster than the temporal evolution of the complex 
valued amplitude A{t) of the oscillations. Therefore the 
oscillation may be written as a product of both factors, 
i.e. oc ^(t)e"^'=*, and there is a very general approach, 
the so-called amplitude expansion, for separating the dy- 
namics at these two disparate time scales [27,28]. The 
amplitude equation describing the evolution of the am- 
plitude A{t) is obtained by a perturbation expansion of 
the reaction equations with respect to the slowly varying 
amplitude A{t) and it is of the form 

TodtA = e{l + ia)A-g{l + ic)\A\^ A. (1) 

The control parameter e measures the relative distance 
from the bifurcation point and tq is the relaxation time 
defined by tq = e/ur-, that depends on the system. If 
the coefficient g of the nonlinear term is positive, the bi- 
furcation to the oscillatory state is supercritical and if it 
is negative, the bifurcation is subcritical. The imaginary 
parts of the prefactors describe the linear and the nonlin- 
ear frequency dispersion. Especially about the extension 
Eq. (1) including spatial degrees of freedom, there exists a 
rich literature as summarized e.g. in a recent review [29]. 
Here in this work we calculate the coefficients of the uni- 
versal equation (1) for microtubule polymerization and 
we discuss their variation in terms of the reaction rates. 



In the following section II we describe the main steps of 
the reaction cycle for microtubule polymerization and the 
respective equations for two models are presented. The 
time-independent solutions for the stationary polymer- 
ization are given for both models analytically in section 
III. Those become unstable against oscillatory perturba- 
tions in the range of high of tubulin-dimer density. The 
respective linear stability analysis and the derivation of 
the oscillation threshold are given in section IV, includ- 
ing their dependence on the reaction parameters. Read- 
ers who are mainly interested in numerical results about 
the oscillation threshold may proceed directly to section 
IV A 4. The partial differential equations for growing and 
shrinking microtubules are of first order in the length mi- 
crotubules and first order in time. Their straight forward 
discretization and numerical solution has to be consid- 
ered with care, therefore a stable numerical scheme, that 
becomes exact close to the threshold of the Hopf bifur- 
cation, is described in section V. The derivation of the 
universal equation (1) is outlined in Sec. VI whereas the 
technical details are given in appendix A. With a sum- 
mary and an outlook about modeling microtubule poly- 
merization we conclude this work in Sec. VII. 



II. MODELS FOR MICROTUBULE 
POLYMERIZATION 

Microtubule assembly and disassembly proceeds in sev- 
eral steps [1-3,5,7,8]. Aggregation of guanosine triphos- 
phate (GTP) liganded tubulin dimers, the so-called 
tubulin-t, to microtubules is started by heating up tubu- 
lin solutions to a temperature of about 30 — 37°C in the 
presence of GTP. Then microtubules spontaneously nu- 
cleate and polymerize to long rigid polymers made up of 
a — (3 tubulin dimers. An increasing number of long mi- 
crotubules in a solvent causes an increasing turbidity and 
the amount of polymerized tubulin-t may be monitored 
by measuring this turbidity [6] or by X-ray scattering 
[8]. The nucleation of microtubules is a rather complex 
process and it is still a matter of debate whether the nu- 
cleation rate depends in experiments only on the initial 
concentration of tubulin-t, Cj, or during the polymeriza- 
tion on the temporally varying Ct [3,30]. But once micro- 
tubules are formed, they grow and the available tubulin- 
t dimers will be used up. The growth velocity of mi- 
crotubules, fg, is rather sensitive to temperature varia- 
tions but it is rather independent on Ct [30,31]. Growing 
microtubules may change their state to rapidly depoly- 
merizing ones by the so called catastrophe rate feat- In 
previous works for the catastrophe rate mostly an expo- 
nential dependence on the tubulin-t concentration was 
assumed, i.e. feat ^ cxp{—ct/cf) with some constant Cf 
[6,15]. Once microtubules have changed from growth to 
shrinking, they shrink rather quickly with a large velocity 

Vs > Vg. 

During the depolymerization of microtubules they are 
fragmented into oligomers or directly into guanosine 
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diphosphate (GDP) ligandcd tubulin dimcrs, the so- 
called tubulin-d dimers. The oligomers themselves are 
believed to fragment further into tubulin d dimers and 
the decay rate depends on the free GTP and GDP. 
Oligomers are stabilized by GDP and destabilized by 
GTP [8,11]. If an excess of GTP is available, then 
tubulin-d in solution will exchange its unit of GDP for 
GTP and each tubulin-t dimer resulting from such an 
exchange step is identical to the initial tubulin-t dimer. 
Such a regeneration step completes the whole micro- 
tubule polymerization cycle. If a continuous source of 
GTP is provided, for instance by a regeneration process, 
this cycling may be continued over a long time [19]. The 
variation of the reaction rates of the polymerization cy- 
cle with the concentration a may depend on the specific 
experiment. 
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FIG. 1. Two models for the cycle of microtubule polymer- 
ization. Model I (upper cycle): Tubulin-t dimers may spon- 
taneously form nuclei of microtubule that grow further by 
incorporating tubulin-t dimers. A growing microtubule may 
also change its state to a quickly depolymerizing one by the 
so called catastrophe rate feat, but it may also change back 
to the polymerizing state by a so-called and rather small res- 
cue rate fresc- Tubulin-d dimers are released during this 
microtubule depolymorization and the whole cycle becomes 
closed by regenerating them by a rate a back to tubulin-t 
dimers. Model II (lower cycle): Here the intermediate step of 
shrinking microtubules is replaced by oligomers e.g. micro- 
tubule break off with a rate feat directly into oligomers and 
the ohgomers themselves may break off with the rate x into 
tubulin-d dimers. The rest of the cycle is identical with the 
upper cycle. 



There are rather detailed models available to describe 
this reaction cycle of microtubule polymerization, see e.g. 
[15]. As a simplification of this complex biochemical re- 
action we only take into account as rate limiting factors 
one of the two intermediate steps of the polymerization 
cycle, either the dynamics of shrinking microtubules (up- 
per cycle in Fig. 1) or the decay dynamics of oligomers 
(lower cycle in Fig. 1). Without both rate rate limit- 
ing factors there are no microtubule oscillations, but one 
is sufficient for oscillations. The two simplified reaction 
schemes, as sketched in Fig. 1, are analyzed in detail in 
this work. 



A. Dynamics of growing microtubules 

During microtubule polymerization there arc many 
growing filaments in a unit volume and their length distri- 
bution may be described by a length and time-dependent 
function Pg{l,t) whose detailed form varies with the ex- 
perimental conditions. A simple model for the dynam- 
ics of distribution of growing microtubules Pg{l,t) is de- 
scribed by the following first order differential equation 
[13,14] 



9tPg = -fcatPg 



dpg_ 
' dl 



(2) 



feat describes cither the transition from the growing to 
the shrinking state of microtubules (model I) or the de- 
cay of growing microtubules into oligomers (model II). 
Vg is the growth velocity of the microtubules. 



1. Growth velocity and catastrophe rate 

In recent experiments with high tubulin-t concentra- 
tion the growth velocity Vg was rather independent of 
Ct [30]. Since we are mainly interested in the oscilla- 
tory behavior of microtubule polymerization, that occurs 
at high Ct - concentrations, we assume a constant Vg in 
this work. In most of the present models a Ct-dependent 
catastrophe rate feat is crucial for oscillatory polymeriza- 
tion of microtubules. Rather common is an exponential 
Cj-dependence [15] 



fcat{Ct) = / e" 



-ct/cf 



(3) 



with the amplitude / and the decay constant c/. How- 
ever, also a linear Ct-dependence 



fcat{Ct) = / (Cu - Ct) ; 



(4) 



with an appropriate constant c„ > Ct leads to an oscilla- 
tory microtubule polymerization as we show in Sec. IV A. 
A hyperbolic Ct -dependence of feat, as discussed in 
Ref. [16], also supports oscillating polymerization. 



3 



2. Nucleation and boundary conditions 

The nucleation process of microtubules is rather com- 
plex and it has been investigated in more detail in Rcfs. 
[12,31,32], recently. The nucleation rate v depends on the 
initial concentration co of tubulin dimers, but it is rather 
independent of the temporal variation of ct, as observed 
in recent experiments [30,3]. Accordingly, for a given 
initial concentration cq we assume a constant nucleation 
rate v. The nucleation rate u itself defines a boundary 
condition for the length distribution of growing micro- 
tubules Pg{l,t) at Z = 0, 



B. Model I includes the dynamics of shrinking 
microtubules 

In model I we take into account as an intermediate step 
between growing microtubules and tubulin d dimers the 
dynamics of shrinking microtubules, Ps{l,t). Here the 
catastrophe rate feat describes the transition of micro- 
tubules from the growing to the shrinking state. The 
depolymerization speed Vs of shrinking microtubules, 
Ps{l,t), is mostly much larger than the growth velocity 
Vg. Having microtubules in two different states, one may 
also expect a transition from the shrinking back to the 
growing state, as described by a rate fresc- So one has 
two coupled equations for the growing and shrinking mi- 
crotubules [13,14] 



dtPg 

dtPs 



-fcatPg + fresc Ps " VgdiPg 
fcatPg - fresc Ps + V^dlPs . 



(6a) 
(6b) 



The rescue rate fresc, however, is usually very small in 
experiments and therefore it is neglected in this work. 
The boundary condition for shrinking microtubules is 



Ps{l ^ oo,t) = 0, 



(7) 



because the transition from growing to shrinking micro- 
tubules is the only source for the shrinking ones and 
Pg{l — > oo, t) vanishes for large values of I. 

The temporal evolution of the concentration of the 
tubulin-t dimers Ct and tubulin-d dimers Cd is described 
by two equations as follows 



+ acd.. 



dtCt = -^Vg I dlPg{l,t) 

Jo 

POO 

dfCd = IVs I dlps{l, t) - acd ■ 
Jo 



(8a) 
(8b) 



The first term in Eq. (8a) describes the consumption of 
tubulin-t during the polymerization (growth) of micro- 
tubules and 7 is a length factor describing the number 



of tubulin dimers that are incorporated in a unit length 
of microtubules. Cj is regenerated from Cd by exchange 
the unit GDP for GTP and this regeneration process, 
described by the rate a, occurs in Eq. (8a) as a source 
and in Eq. (8b) as a sink. Tubulin-d dimers are released 
during the depolymerization of microtubules Ps{l,t) and 
this source is described by the integral in Eq. (8b). 

Tubulin dimers may be a constituent of growing or 
shrinking microtubules or they carry GTP or GDP as 
single dimers, but altogether they are conserved as ex- 
pressed by the condition 



ct + Cd + 'yL = Co. 



(9) 



Here cq describes the overall concentration of tubulin 
dimers and L{t) is the integrated length of all micro- 
tubules per unit volume 



L{t) = J^ dl l{pg{l,t)+ps{l,t)y 



(10) 



The tubulin d concentration may be eliminated from 
Eq. (8a) by using the conservation condition (9). On 
the other hand Eq. (9) in combination with Eqs. (6) and 
Eq. (8a) yield an equation that is identical with Eq. (8b). 
Hence Eqs. (6a) and (6b) together with 

dtCt ^ ~'y j '^^ i^gPg + ^KPg +Ps)) + "(co - c*) , (11) 



describe the polymerization dynamics of microtubules for 
model I, whereby a constant growth and shrinking veloc- 
ity is assumed in this work. 



1. Resettling of model I 



After rescaling time t and length I, i.e. 



t' = at, 



(12) 



it is easy to see that model I may be characterized by a 
set of dimensionless parameters 



a 



a 



Co f, 

Cf 



a 



(13) 



Some of these dimensionless quantities may be further 
combined to other dimensionless parameters as for in- 
stance in the threshold condition given in Sec. IV A. 
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2. Reduced model 

Since the depolymerization velocity Vg is much larger 

than the growth velocity Vg one may also consider the 
limit Vs ^ Vg. In this case the shrinking micro- 
tubules decompose nearly instantaneously into tubulin-d 
dimers and growing microtubules decay effectively, due 
to the short life time of the shrinking microtubules, into 
tubulin-d dimers. In order to describe this direct decay 
the source term in Eq. (8b), 71^5 /q dlps{l,t), must be 
replaced by 7/ca* Jq^ dllpg{l,t). Eliminating again the 
density Cd one ends up with a reduced model for only 
two densities: 

dtPg = -fcatPg - VgdlPg , (14a) 

dtCt = ~T y '^^ (^9 + "^)ps + Q'(co - ct) ■ (14b) 

This simplified model reproduces essential aspects of sta- 
tionary polymerization of microtubules as described in 
Sec. III. 



C. Model II includes the dynamics of oligomers 

Oligomers occur as an intermediate product during the 
decay of microtubules and they are made of several tubu- 
lin dimers. This intermediate product is ignored in model 
I. Here in model II, after the so-called catastrophe, we 
ignore the dynamics of shrinking microtubules as an in- 
termediate step and instead we take into account the (de- 
cay) dynamics of oligomers. Therefore, the catastrophe 
rate feat in Eq. (2) describes for model II a direct tran- 
sition of growing microtubules into oligomers. Further- 
more, it is assumed that oligomers decay with the rate x 
into tubulin-d dimers. The concentration of oligomers is 
denoted by Cqh and its dynamics as well as that of q are 
described by the two equations 

POO 

dtColi = rjfcat / dllpg{l, t) - xcoli , (15a) 
Jo 

dtCd = xA Coii - acd ■ (15b) 

ry is a measure for the number of oligomers per unit length 
of the microtubules and A is a measure for the number 
of tubulin dimers per oligomer. Oligomers decaying with 
the rate x build a source term in the equation for tubulin- 
d dimers in Eq. (15b). 

The conservation law for the concentration of all tubu- 
lin dimers takes the form 



Eq. (2) and Eq. (8a) together with the following dynam- 
ical equation for Cd 



r 

Ct+Cd + XcoH +ri\ I dll pg{l, t) = Co. 
Jo 



(16) 



The equation for the growing microtubules is the same 
as for model I, cf. Eq. (2), but in the equation for Cj, 
cf. Eq. (8a), one has to replace the length factor 7 by 
the product rjX. Eliminating Cqh model II is described by 



dtCd = x[co-Ct-Cd-v>' J dl Ipgj - acd ■ (17) 

As boundary condition for the growing microtubules we 
again use Eq. (5) with a constant nucleation rate v. For 
model II wc only consider the catastrophe rate given in 
Eq. (3). This again guarantees a nonlinear feedback of 
the dynamics of the tubulin-t dimers to the dynamics of 
the growing microtubules. 



1. Reduced model 

Similar as for model I also model II becomes in the 
limit X — > cx) identical with the model described by 
Eqs. (14). If we assume a very fast dissociation of 
oligomers into tubulin-d dimers, x ^ 1) we can neglect 
the intermediate state CoU- In this case the source term 
in equation (15b), xAcofi, can be replaced by the source 
in equation (15a), cf. rjXfcat /o°° dllpg, which describes 
the direct decay of growing microtubules into tubulin-d 
dimers. After replacing Cd and setting 7 = rjX in Eq. (8a), 
we again obtain with the help of the; conservation law (16) 
the simple reduced model as described by the equations 
(14). 



III. STATIONARY SOLUTIONS 

A polymerization cycle with a stationary length distri- 
bution of microtubules and time-independent dimer con- 
centrations Ct, Cd or oligomer concentration Coh are one 
type of the solutions of the model equations described 
in the previous section II. For this stationary state the 
various polymerization steps, such as nucleation, assem- 
bly and disassembly of microtubules as well as the re- 
generation of tubulin-d dimers are in a balanced state. 
Under certain conditions a stationary polymerization is 
observed in experiments [31]. However, it may become 
unstable against oscillatory perturbations if the initial 
tubulin dimer concentration cq is large enough, as shown 
in the following section IV. 



A. Model I 

Eqs. (6) are first order linear differential equations with 
respect to the length I and in the stationary case these 
equations have exponentially decaying solutions, which 
take in the limit of a vanishing rescue rate the form 



^>^?l(0 = -^expf-A 



(18) 
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,(0) 



t concentration, denoted by c, 
consistently as described below. 

JO) 



is determined self- 

The stationary solu- 



The catastrophe rate /^°| may be given either by Eq. (3) microtubules which consist of many tubulin dimers and 
or Eq. (4) but in both cases the stationary tubulin- therefore enhance the densities ^ and c^°^ . The length 

distribution of the growing respective shrinking micro- 
tubules is determined by Vg/ f'^^, which also depends via 
the catastrophe rate on the concentration Cj""* . 

These tendencies becomes even more obvious if the lin- 
ear dependence of the catastrophe rate in Eq. (4) is cho- 
sen for the special case c„ = co and Eq. (19) is expanded 
in the limit of small and large values of a. In both cases 
we obtain the simple formulas 



tions pg.s allow an analytical calculation of the integrals 
in Eq. (11) and a nonlinear equation in cf^ follows, 

:(1 + /3)1 • 



Co 



AO) 

J cat 



1 

7(0)' 

./ cat 



(19) 



From this equation the stationary tubulin concentration 
c[^^ can be determined as a function of the overall con- 
centration of tubulin dimers cq and as function of the 
other parameters. In Eq. (19) the abbreviation for the 
velocity ratio 



(20) 



has been introduced and the respective length distribu- 
tions pf ^ and p^'s'^ follow for a given value of c^^^ via 

Eqs. (18). The stationary value cf^^ for the reduced 
model, described by Eqs. (14), follows from Eq. (19) in 
the limit /? ^ 0. 
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FIG. 2. The tubulin-t concentration for the stationary 

polymerization state of model I is shown as a function of the 
regeneration rate a and for two different values of the nu- 
cleation rate v. The velocity ratio between the growing and 
shrinking microtubules is /3 = Vgjvs =0.1 and the rest of 
parameters are co = 120, Vg = 0.1, c/ = 3, / = 0.1, 7=1. 
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which reflect the described tendencies. 



B. Model II 



Stationary solutions for model II can be calculated in 
a similar manner as discussed in the previous section for 
model I. The length distribution of growing microtubules 
is again given by Eq. (18) and the integral in Eq. (17) 
can be calculated analytically, cj. may be eliminated 
from Eq. (17) by using Eq. (8a) with ct = and by set- 
ting 7 = 77A. Then the nonlinear equation for tubulin-t 
(0) 



dimers c^^' takes the form 



Co 



,(0) 



vrjXvg 

J cat 



Eq. (22) is invariant under permutation a ^ x- If o; and 
X become much larger than the catastrophe rate, the sta- 



becomes rather independent of 

00 in Eq. (19) and x ^ 00 in 

(0) 



tionary concentration Cj 
both. In the limits v, 

Eq. (22) we again obtain the concentration c:^'' for the 
reduced model. No stationary solution is possible in the 
limit X ^ because in this limit all tubulin-d dimers are 
stored in oligomers and the polymerization cycle becomes 
interrupted. 



In the range of a much larger than the catastrophe rate 
/^^°| the stationary tubulin-t concentration c|"-' becomes 
independent of it, because all tubulin d dimers, that are 
released during the depolymerization, are immediately 
regenerated to tubulin-t dimers. Both a large nucleation 
rate v and a large growth velocity Vg lead to a high con- 
sumption of tubulin-t and therefore to a lower station- 
ary concentration cf'\ This tendency is illustrated by 
the difference between the two curves in Figure 2. On 
the other hand, large values of the amplitude of the re- 
spective catastrophe rate, either / or /, act against long 



IV. THRESHOLD FOR OSCILLATORY 
POLYMERIZATION 

Stationary microtubule polymerization becomes unsta- 
ble against oscillating modes in the range of high tubulin 
dimer concentrations cq and the parameter range where 
this happens is calculated by a linear stability analysis. 
Starting from the model equations given in Sec. II we de- 
rive linear equations for small perturbations with respect 
to the stationary state and such perturbations exhibit 
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an exponential time dependence, ef^. For the exponen- 
tial factor a we derive a nonlinear equation from which 
both the critical dimcr concentration cqc and the criti- 
cal frequency lOc for the Hopf bifurcation is calculated 
numerically for various parameter combinations and in 
limiting cases also analytically. 



Herein we have introduced the abbreviations 



(0) 



(0) ' 



" J cat 



(28) 



A. Model I 



We introduce small perturbations p'^gj and c[^^ with 
respect to the stationary solutions pf} and c[^^ as deter- 



mined in the last section. With the ansatz 
Pg,s 



Pg,s > 



„ _J0) , „(1) 
Ct — C( -|- Cj 



(23a) 
(23b) 



one obtains, after linearization of Eqs. (6) and Eq. (11), 
the following set of linear equations describing the dy- 
namics of the perturbations 



dtp^^ = -{f!^+vgd;)pi^^-p^^f^, 



(24a) 
(24b) 



(1) 



„(i) 



. dl + a l{p^^^ + p(i))) . (24c 



Here /^^^ is the first order contribution of an expansion of 



the catastrophe rate /c, 



to the perturbation q 



(1). 



AO) 

J cat 



J cat 



with respect 



and the boundary condition in Eq. (7) requires a vanish- 
ing integration constant K = 0. The boundary condition 
for the time-dependent part of the growing microtubules, 
p^g\l = 0,t) = 0, is also fulfilled. According to the ana- 
lytic expressions for pg^-* and pi^-* given in Eqs. (27) both 
may be eliminated in Eq. (24c). The remaining integral 
in Eq. (24c) can be calculated analytically and the non- 
linear dispersion relation for a follows 



(29) 



AO) 

J cat 



■f3- 



AO) 

J cat 



/i°)+a(l+/3)^ 



with a reduced parameter 



G = 



(30) 



for the catastrophe rate given in Eq. (3) and with 

AO) 

G= 4^ (31) 



Ai) 

J cat J ( 



(0) 



(25) 



Since the first order linear equations (24) have constant 
coefficients, their solutions depend exponentially in time 
and C(^^ may be written as 



=Ae''* + c.c. 



(26) 



(c.c.=conjugate complex). With this ansatz the three 
equations in (24) can easily be integrated and the so- 
lutions for the growing and shrinking microtubules are 
given by 



.(1) 



K[cat_ 

VgCfa 



exp I at 



exp 



(0) 

cat 



VgCfG 



exp at 



AO) 

J cat 



A 



AO) 

J cat 



C.C. 



(27a) 



+ k2 + K ■ exp 




(27b) 



for the rate given in Eq. (4). After a few rearrangements 
the dispersion relation in Eq. (29) can be written as a 
fourth order polynomial in a 



a^ G0 {l+p)+a^G ap{l + (3) + f^°l {1 + 3(3 + (3^) 



+a^ 



Gaf!:°l (1 + 3/3 + (3^) + (/? + 2Gf!,f) (1 + (3) 
a{l + f3) (l + /3 + 2G/i°r)+G/, 



AO)-" 

I cat 



+ (1 + 2^) 



2 + 2/3 + G/i°r 



(32) 



This polynomial describes the linear stability of the sta- 
tionary solutions given by Eq. (18) and Eq. (19) com- 
pletely and they are unstable in the parameter range 
where the growth rate becomes positive, Re{a) > 0. 
Keeping for instance all parameters besides the dimer 
concentration cq fixed, then the neutral stability condi- 
tion Re{a) = provides an equation for the critical dimer 
concentration cqc- For concentrations larger than this 
critical value, cq > cqc, the stationary solutions are un- 
stable. 

The smallest critical dimer concentrations cqc for an 
oscillatory polymerization are required if the parameters 
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a,/3 and G take intermediate values, as discussed in more 
detail below. At the threshold the real part Re{(T) = 
vanishes and the imaginary part of a is the so-called 
Hopf-frequency lOc = Im{a). In this special case with 
a purely imaginary a the polynomial in Eq. (32) can be 
decomposed into its real and imaginary parts giving two 
coupled equations for the determination of the two un- 
knowns and ujc- Having determined /^'^^ numerically, 
coc may be calculated via Eq. (19). 



2. Limiting cases for the rate in Eq. (4) 

The tendencies for the parameter dependence of 
the threshold for the Hopf bifurcation as discussed in 
Sec. IV A 1 are by far not a special property of the choice 
of the catastrophe rate in Eq. (3). Therefore we consider 
the same limiting cases as before for the catastrophe rate 
given in Eq. (4) and with G as defined in Eq. (31). 

a. a — > oo; In this limit one obtains from Eq. (32) 



1. Limiting cases with the rate in Eq. (3) 



For the limiting cases /3^0,/3^oo,Q!^0 and 
a — > oo analytical expressions can be given for both the 
threshold concentration cg^ and the Hopf frequency uJc- 
This is explained at first for the catastrophe rate given 
by Eq. (3) and for the parameter G given in Eq. (30). 
At threshold one has a = iu)c and two equations follow 
from the nonlinear dispersion relation in Eq. (32) which 

determine the two unknowns uJc and f^^l- The critical 
initial concentration cqc follows via /^°] from Eq. (19). 
a. a — > oo; In this limit one obtains from (32) 



J cat 



1 + 



aG 1 + /3 + /?2 ' 




(33a) 
(33b) 



AO) 

•I cat 



5(1 + 



1/2 



a(l+/3 + /?2)^ 
(ag(l+/3 + /3^)(l + /?))^/' 



(35a) 
(35b) 



with g = f^vvg. Hence the critical tubulin concentration 
required for a Hopf bifurcation diverges as cqc oc a. 

b. a 0; In this limit one has fj!fl oc a and cqc oc a~'^ 
diverges too. 

c. P ^ 0: For this limit we obtain 



AO) 

J cat 



«l/6„l/2 I ± 



1\ 1/3 



-0) 



(36a) 
(36b) 



This confirms the importance of a finite ratio oi (3 = 
Vg/vs, because the threshold diverges for /? ^ 0, similar 
as for the catastrophe rate given in Eq. (3). 



Accordingly the critical tubulin concentration diverges 
as coc oc a^, that agrees with the full numerical results 
shown in Fig. 3, besides small logarithmic corrections. In 
this limit the Hopf frequency lOc becomes independent of 
a and with increasing values of /? it decreases slightly to 
a constant value Uc ^ \f\jG. 

b. a ^ 0.' In this case lo,, ^ yl/G becomes also 
independent of a and the catastrophe rate vanishes as 
feat ^- Therefore the critical tubulin concentration 
diverges according to Eq. (19) as cqc ~ a~^. 

c. (3 0: In this limit one obtains 



AO) 

J cat 



1/4 



1/4 



(34) 



The Hopf frequency diverges with increasing values of 

the shrinking velocity as ujc ^ vl^^ in agreement with 
the numerical results shown in Fig. 4. With this ex- 
pression for one obtains via Eq. (19) for the critical 
initial concentration cqc ^ G//? + cy In (/(G//3)"^/^). For 
medium parameter values this is essentially cqc oc 1//3 as 
indicated in Fig. 4. 

In experiments the shrinking velocity was always larger 
then the growth velocity, therefore the limit /3 ^ oo is 
discarded. 



3. Traveling waves solutions 

At the threshold of the Hopf bifurcation the rate a is 
purely imaginary, a = iuic, and the expressions given in 
Eq. (26) and Eqs. (27) are oscillatory in time 



,(1) 



2 A cos{wct) , 



rf) = ^exp(-:^0 



(0) 



(37a) 

sin{ujct) — sm{ujc{t — l/vg)) , (37b) 



Si 



AO) 

J cal 



exp(-^-^/) [k2 sin(a;ct + <^2) 



+ ki sin(a;c(i - l/vg) + (fi)] , (37c) 
whereby the following abbreviations for the amplitudes 



Si 
ki 

k2 = 



OJc Cf 

V/iar+^g/if(i+/jF 

/if +C.2(l + /^)2 



('^c'/?-/if) +/i°ra;2(l + /3)^ 
/if + (^c/3)' 



(38a) 
(38b) 

(38c) 
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and phases 



ifii = — arctan 



V?2 = arctan 



^ / cat / 

fi^i - /IT . 



(39a) 
(39b) 



Hopf bifurcation for the reduced model that follows in 
the limit /3 — ^ as described in section II B 2. Hence, 
the dynamics of shrinking microtubules is one essential 
degree of freedom favoring oscillating microtubule poly- 
merization. The dynamics of oligomers, as discussed in 
Sec. IV B, is an alternative degree of freedom that favors 
oscillations. 



have been introduced. The analytical expressions for p^p 
and p^p indicate that the time-dependent contribution 
to the length distribution of the microtubules includes 
homogeneous amplitude oscillations and waves with a 
wavelength ojc/vg that travel to larger values of the length 

I. Hence, the length distributions p^g}s depend on two 
different length scales, the decay length Vg/f^^l and the 
wave length Vg/cUc of the traveling waves. With the ex- 
plicit solutions for c[^^ and Pg^' the phase of the oscil- 
lating part of the tubulin-d concentration relative to the 
phase of Cj^^ as well as its oscillation amplitude is calcu- 
lated via Eq. (8a). 



4- Numerical results for the threshold of model I 

Since a is purely imaginary at the threshold of a Hopf 
bifurcation, a = itOc, the dispersion relation in Eq. (32) 
can be decomposed in its real and imaginary part and 
from these two equations the critical concentration coc 
and the frequency Uc may be calculated numerically as a 
function of the parameters. The numerical calculations 
in this section are restricted to the exponential tubulin 
dependence of the catastrophe rate as given by Eq. (3). 

The critical tubulin dimer concentration cqc and the 
critical frequency ojc at the Hopf bifurcation are shown in 
Fig. 3 as function of the regeneration rate a and in Fig. 4 
as function of the velocity ratio = Vs/vg, whereby 
the reduced parameter G has been chosen at the values 
G = 3000 and G = 300, respectively. Since G includes 
a number of parameters the curves in both figures repre- 
sent a larger parameter set. In the limit of a vanishing 
regeneration and in the limit of very large values of a, 
where the regeneration process is much faster than any 
other process, the critical tubulin concentration cqc di- 
verges and therefore the Hopf bifurcation is suppressed. 
In addition, both figures indicate that the smallest values 
of the critical tubulin concentration cqc are obtained at 
intermediate values of the parameters a, j3 and G. The 
location of the threshold minima, however, depends on 
the actual values of the rest of the parameters. The fre- 
quency ijjc becomes rather small in the limit a ^ and 
for large values of a this frequency becomes independent 
of it, cf. section IV A 1. 

In Fig. 4 the threshold minimum is less pronounced 
than in Fig. 3 and in the limit (3 = Vg/vs the thresh- 
old coc increases linearly with Vg and in agreement with 
limits given in section IV A 1. Accordingly, there is no 



0.06 



o 

(D 



g 




O 

c 
o 
o 




0.15 



regeneration rate a 

FIG. 3. The critical tubulin dimer concentration coc and 
the critical oscillation frequency oJc are given at the threshold 
of the Hopf bifurcation as a function of the regeneration rate 

a, for two values of /3 = Vg/vs and for the constant G = 3000. 
The catastrophe rate given in Eq. (3) has been used with the 
parameter values / = 0.1 and c/ = 3. 

For large values of Vg the frequency cjc becomes large 
too and the oscillation period becomes much shorter than 
any rclaxational dynaniic;s of and Cf According to the 
quick shrinking, the life time of a depolymerizing micro- 
tubules vanishes and therefore the amplitude of the den- 
sity of shrinking microtubules is small too, ps ac l/vs- In 
other words, in the limit of large values of Vs the inter- 
mediate step of shrinking microtubules may be neglected 
and the transition from pg to tubulin-d dimers is effec- 
tively a direct process as explicitly assumed for the re- 
duced model. If either the regeneration or the shrinking 
dynamics becomes too fast, the Hopf bifurcation is sup- 
pressed. The two intermediate steps, the depolymeriza- 
tion and the regeneration, act obviously as antagonistic 
steps or jam processes that favors oscillations. 

Since one has at threshold a = iu>c, Pg^^ and pi^^ 
in Eq. (37b) and Eq. (37c) include both traveling wave 
contributions cx exjp[—i(LUct — kl)] with a wave number 
k = uJcjVg and that always travel towards larger lengths 
of the microtubules. The length distribution is exponen- 
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(0) 



If this is large 
as for instance 



tially decaying on the length scale Vg / /, 
compared to the wave length A = 2'KVg/u 
in the limit Vg Vg, then one has a kind of self averag- 
ing in the respective integrals and the Hopf bifurcation 
is suppressed. 
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velocity ratio P = v / v 



g 



FIG. 4. The critical tubulin concentration coc and the fre- 
quency ujc are given at the Hopf bifurcation as a function of 
the ratio between the shrinking and growth velocity of mi- 
crotubules, Ve/vg = (3~^ , for a = 0.05 and for two different 
values of G. The ct-dependence of the catastrophe rate as 
given in Eq. (3) has been used with the same parameters as 
in Fig. 3. 

The phase difference between the oscillations of 
tubulin-t and the oscillations of the total amount of poly- 
merized tubulin, described by L{t) in Eq. (10), is another 
experimentally accessible quantity [33]. The difference 
between the phases of the oscillatory contributions of 
and Ct as well as the difference between the phases of 
L{t) and Ct are given in Fig. 5. These phase differences 
as well as the ratios between the amplitudes of the fields, 
cf. lower part of Fig. 5, are calculated at the threshold 
of the Hopf bifurcation by using the analytical solutions 
calculated in Sec. IV A 3. 

For large values of a, tubulin-d is quickly regener- 
ated into tubulin-t and therefore the density q becomes 
smaller as shown by the lower part in Fig. 5. In the op- 
posite limit of small values of a tubulin t is consumed by 
nucleation and growth of microtubules, but the source, 
which is supplied by the regeneration of tubulin-d, de- 
cays and therefore one obtains large values for the ratio 



between the amplitude of c^^'' and c^^' as well as between 
the amplitudes of i^^^ and c^^\ 

The decay of the ratio between the amplitudes of L^^^ 



(1) 



and is less obvious, 
contributions L^g^ 



= 7/0" 
(1) _ 



dll 



Acos{uJct~\'(p^) and 
L^p = 7 dllpP = Bcos{aJct + ipg). The two am- 
plitudes A and B increase with the regeneration rate a. 
However, with increasing values of a the phase difference 
(fix — ip§ increases as well up to unity leading to an ef- 
fective decay of the sum L^^^ as shown by the lower part 
of Fig. 5. 
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FIG. 5. In the upper part the differences between the 
phases of the oscillating contributions of the tubulin-d con- 
centration Cd (solid) and of the tubulin-t Ct as well as between 
L{t) and ct (dashed) are shown as function of the regeneration 
rate a. In the lower part the ratios between the amplitudes 



of the oscillating contributions and cj^' (solid) as well as 
between the amplitudes of L*-^' and c[^^ (dashed) are shown. 
The rest of parameters are as in Fig. 3. 



The phase shifts of 04 and L{t) with respect to the 

phase of Cj are rather independent of the regeneration 
rate a as shown in Fig. 5. The absolute values of these 
shifts are in qualitative agreement with the expectation 
as described in the following. At the maximum of ct 
the catastrophe rate takes its minimum and therefore, 
since the nucleation and the growth velocity are constant, 
L{t) is increasing for a while and up to the moment when 
enough Ct is consumed and the catastrophe rate increases 
again. Due to an increasing decay of microtubules the 
maximum of the latter will also lead to a delayed max- 
imum for Cd- For large values of a the amount of poly- 
merized tubulin is nearly in anti phase with respect to Cf , 
which is also mentioned in Ref. [15]. The slightly stronger 
a-dependence of the phase difference between L{t) and 
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Cf is mainly due to the a-dependence of shrinking micro- 
tubules Ps, because the relative phase of p^P is nearly 
independent of a (see also Section IV B). 



B. Model II 

Hero the stability of the stationary polymerization 
state of model II, described by cf'\cf^ and pf\^ is in- 
vestigated with respect to small perturbations c'^\ c^P 



and p^g\ With the ansatz 



(0) 



(40a) 
(40b) 



the equations for model II are linearized with respect to 
these perturbations and one obtains the following set of 
linear differential equations with constant coefficients 



poo 



.(1) , .(1) 



/g^] is the first order correction with respect to its value 
in the stationary state and it is given by equation (25). 

The time -dependent contributions to the tubulin-t 
and tubulin-d dimer densities are described by 



ac 



(41a) 
(41b) 

(41c) 



=Ae''' + c.c. , 



EAe''* +C.C. , 



(42a) 
(42b) 



with the common complex amplitude A and relative com- 
plex factor E that describes via E = \E\e^'^'^ the ampli- 
tude ratio \E\ as well as the phase difference ipd between 
both fields. With the solution for growing microtubules 



as given in Eq. (27a) we can eliminate c^^^ and Pg"^' in 
Eq. (41b) and we again obtain from the resulting sol- 
ubility condition a nonlinear dispersion relation for the 
exponential factor a 



+ ax 




(43) 



with G = Cf / {rjXvgv). After a few rearrangements of this 
equation one obtains a fourth order polynomial in a for 
model II as well 



-Vi°lG + a3/i°)G(2/i°)+a + x) 
+ [l + ocxG + G/i"] (2x + /iS + 2a) 

+ (/i°i+« + x)+aX 

+ G(2ax + /i°|(a + x))" 

+ /iS ax {Gf^f + 2)+ /if (a + x) = , (44) 

that determines the linear stability of the stationary poly- 
merization for model II. Again we are interested in the 
neutrally stable case, Re{a) = 0, that separates the sta- 
ble from the unstable regime. At the neutral stability 
point of the Hopf bifurcation one has Wc = Im{a) and 
Eq. (44) can be decomposed into its real and imaginary 
part. From these two equations and ujc are deter- 
mined by standard methods, cqc may be calculated via 
Eq. (22). 



1. Traveling waves solutions 

At the Hopf bifurcation the non-stationary part of 
growing microtubules is again described by the distribu- 
tion given by Eq. (37b) and the fields Cqu and are not 
in phase with cj in general. The two fields may be writ- 
ten in terms of the amplitude ratio |£^| and the relative 
phase ^pd in the following form 



^ = 2 A cos(wci) , 



(45a) 
(45b) 



The amplitude ratio \E\ and the phase ip^ can be deter- 
mined from the two coupled equations (41b) and (41c) 
and they are given by 



\E\ = 



ifd = arctan 



\ J cat 




(46) 



In a similar manner the oligomer density c^|] may be 
written in terms of an amplitude ratio \F\ and a relative 

phase ipou between c^)- and Cj^^ 



4)] =2A\F\ cos(wci + ^Poii) , 



(47) 



with 



\F\ 



X 

arctan 



E\ , and 

atan(yrf) -|- Uc 

a - LUcta.n{ipd) 



(48) 



The oscillatory contribution to the polymerized tubu- 
lin can also be written as a harmonic function, 
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L(i) = 2A\H\cos{ i^ct + ^h)- For both, the amplitude 
ratio \H\ and the relative phase </5l, one obtains long ex- 
pressions that are not presented here. The phase shifts 
and the amplitude ratios between the oscillating fields 
are shown in Figure 7 as function of the regeneration 
rate a. As discussed in Section IV A, the phase shift of 
the polymerized tubulin L'^^^(t) with respect to ci^ is 
rather independent of a, whereas the phase of oligomer 
oscillations change slightly with a. A phase shift tt be- 
tween the polymerized tubulin and oligomers is measured 
in experiments, cf. Refs. [12] and [34]. In this model this 
is only possible in the limit of a dissociation rate x much 
smaller than the regeneration rate a. 



2. Numerical results for the threshold of model II 

At the threshold one has again a = icOc and from the 
imaginary together with the real part of the dispersion 
relation in Eq. (44) the critical concentration cqc and the 
Hopf frequency uic niay be calculated as function of the 
parameters. Also for model II we restrict our numerical 
calculations to the catastrophe rate with the exponential 
dependence given in Eq. (3). 
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FIG. 6. For model II the critical tubulin concentration coc 
and the critical frequency ijJc are shown at the Hopf bifur- 
cation as a function of the regeneration rate a. and for two 
different values of the decay rate x of oligomers. The rest of 
the parameters are G = 3000, / = 0.1 and c/ = 3. 

The critical tubulin concentration cqc and the critical 

frequency Lu'c at the Hopf bifurcation are shown in Fig. 6 
as function of the regeneration rate a and for two differ- 
ent values of the decay rate of oligomers Xj whereby for 



the reduced parameter G the value G = 3000 has been 
chosen. Since G includes a number of parameters the 
curves in both parts represent a larger parameter set. 
For a fixed finite value for x in the limit of a vanish- 
ing regeneration a — > 0, where the polymerization cycle 
is interrupted, and in the limit of very large values of 
a, where the regeneration process is much faster than 
any other process, the critical tubulin concentration cqc 
diverges similar as for model I and therefore the Hopf 
bifurcation is suppressed. If a is kept fixed at a medium 
value the threshold curve coc(x) as function of the decay 
rate x for oligomers has a similar shape as shown as func- 
tion of a in Fig. 6. The critical tubulin concentration cqc 
also takes its smallest values at intermediate values of a, 
X and G, whereby the location of the threshold minima 
depends on the actual values of the rest of parameters. 
With a decreasing rate a — > of tubulin regeneration 
also the frequency Wc becomes small. On the other hand 
for large values of a the tubulin regeneration is not any- 
more a rate limiting factor and the critical frequency u)c 
becomes rather independent of a , cf. section IV A 1. 
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FIG. 7. The phase differences (upper part) and the am- 
plitude ratios (lower part) between the oscillating contribu- 
tions to the tubuUn-d concentration c^^' (solid), the oligomer 
concentration c[,)] (dotted) and the total polymerized tubulin 
i'^'(t) (dashed) with respect to the tubulin-t concentration 
Cj^' is shown as a function of the regeneration rate a. The 
other parameters are G = 3000 and x = 0.02. 



The stability of oligomers and therefore the decay rate 
X depend very much on the available GTP: Increasing 
GTP concentrations destabilize oligomers and increase 
the decay rate x [7,11]- With increasing GTP concen- 
trations also the rate a of the transition from ca to C( is 
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enhanced. However, if the tubulin regeneration and the 
oligomer decay become to quick, an oscillatory polymer- 
ization is suppressed. In other words, if one increases a 
and X beyond some minimum values, the threshold con- 
centration for tubulin increases too. Such a tendency for 
the GTP dependence of the oscillation is in agreement 
with the results reported from experiments [7,11,12]. 

With increasing values of a tubulin-d is again quickly 
transfered by the regeneration process into tubidin t, 
leading to a small amplitude ratio c'j^^ / c'l^K Accord- 
ingly more tubulin is left to be stored in L'--^^ and c^j]. 
Therefore both increase with larger values of a as indi- 
cated in Fig. 7. This has to be compared with L^^^ for 
model I, where it decays as function of a because the 
phase shift between the contributions of the growing and 
shrinking microtubules changes too. For model II the rel- 
ative phases are also rather independent of a, whereby 
due to the quick regeneration of Cd the relative phase 
between c^^^ and c^^^ is slightly decreasing. 



where the first mode describes just the stationary so- 
lution and the second one the oscillatory contribution. 
This approximation becomes exact close to the thresh- 
old and this ansatz leads with Eqs. (6) to two differential 
equations for Fg^g 

9tFg = - /cat) + F,) - VgdiFg , (51a) 



+VsdiFs . 



(51b) 



Both fields may be expanded with respect to the first two 
spatial Fourier modes e*""^ (n = 0, 1) 

Fg{l,t) = B{t) + \ [C{tY^^ + C*{t)e-'^^) , (52a) 



F,{l,t) = D{t) + - (i?(t)e*" + i?*(t)e-*") 



(52b) 



3. Reduced models 

The dispersion relation for model I and II, considered 
in the previous section, became equivalent in the limits 
/? ^ and X — > oo and in both cases one obtains the 
same dispersion relation 



cat + 

1 + 2G/,„, 



^ J cat 



+ f 

cat ' J cat 

0, 



(49) 



with the reduced parameter G as given in Eq. (30) for 
model I and with G = Cf/ {rjXvgi^) for model II. This poly- 
nomial in a has always negative growth rates, Re{a) < 0, 
and therefore stationary solutions are always stable. 



V. NUMERICAL METHOD FOR MODEL I 



in order to remove the spatial dependence from Eqs. (51). 
Herein the wave number is chosen at its value at the 
threshold of the Hopf bifurcation, k = iVc/vg. This ansatz 
leads to a set of ordinary differential equations for the 
time dependent amplitudes B{t), C{t), D{t), H{t) that 
are described in the following. 

Due to Eq. (5) one has the boundary condition Fg {I = 
0,t) = that gives the relation Cr = —B, with Re{C) = 
Cfi, between these two functions. Ansatz (52a) together 
with equation (51a) leads to the relation 



Im{C) = Cj = £, {feat - 



(53) 



and to the first order differential equation in B 

dtB^[/f2-feat) (^ + ^) • (54) 

Ansatz (52b) in equation (51b) gives the set of coupled 
differential equations 



The two differential crjuations for growing and shrink- 
ing microtubules in (6) are of first order with respect 
to the length / of the microtubules and first order in 
time. A straight forward spatial discretization of such 
first order equations often leads to numerical instabili- 
ties. Especially the equation for shrinking microtubules, 
cf. Eq. (6b), has problematic stability properties. For 
this reason we approximate the solutions of Eqs. (6) by 
a two mode ansatz 



Pg{l,t) = e^p\-i^l\ \^-+Fg{l,t)] , (50a) 




p,{l,t) = exp\-i^l] i-+F,{l,t) \ , (50b) 



f(0) 



Vs 



Vg J Vg Vg 



cat _ r{0) 

_ . J cat ^ J 



(55a) 



dtHR = -f,at B - ^/(°^ Hr - kVsHi , (55b) 
dtHi = feat Ci - Hi + kvsHR . (55c) 

With the periodic Z-dependence given in Eqs. (52) the 
integrals in Eq. (24c) can be evaluated and one obtains 
the following differential equation for the tubulin-t dimer 
density 

dtct = - jVgK{t) - a'yL{t) + acoc {l + e) - act , (56) 
where the coefficients are given by 
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K{t)= / dlpgil,t) = — 



k{kB - SCi) 



(57) 



L{t) = dll [Pg{l,t)+Ps{l,t)] 

Jo 

V 

+ S^{d^ - k^)HR - 2 kHi + + k'^fD) , (58) 



^ + — ) +A{3d'^k'^B + k^B-25^kCi 

'g 



and where the abbreviations A = [^^((^^ + fc^)^] ^ and 
5 = fc^t/Vg have been introduced. The reduced con- 
trol parameter £ = (cq — cqc) / coc measures the difference 
between the tubuUn dimer concentration cq and the crit- 
ical one, Cqc- For e > sustained oscillations occur but 
they are damped below threshold e < 0. For the numer- 
ical solution of model II we use either the same approx- 
imation scheme, where only the factors of the scheme 
take a different form, or in the absence of Ps a direct 
spatial discretization provides also a stable algorithm. 
The five differential equations for model I in Eq. (54), 
Eq. (55) and Eq. (56) and the corresponding three differ- 
ential equations for model II are integrated numerically 
by a second-order Runge-Kutta method with a time step 
At = 0.01. 
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FIG. 8. For model I the time-dependence of the tubulin-t 
concentration ct{t), the polymerized tubulin jL{t) and the 

tubulin-d concentration Cd{t) is shown in parts a), b) and c), 
respectively. The parameters v — 0.01, a = 0.01, /3 = 0.1 
were used with the corresponding critical initial concentra- 
tion coc = 80.69. The reduced control parameter is chosen 
at the value e = 0.01. In part d) the length distribution of 
the microtubule P(l) = Pg{l) + ps{l) is shown at two different 
times t = 876 (solid) and t = 975 (dashed), where ^L{t) takes 
its maximum and minimum, respectively. 

The time-dependence of the fields of model I are shown 
in Figure 8 for one parameter set and these fields obvi- 



ously have different phases. The extrema (maxima) of 
the polymerized tubulin L(t) and the tubulin d concen- 
tration Cd(t) are delayed with respect to the extrema of 
the tubulin-t concentration Ct{t), a behavior that is al- 
ready indicated by the reaction cycle shown in Fig. 1. At 
the threshold the parameter dependence of these phase 
differences may be calculated from the formulas given in 
Sec. IV A 3 and Sec. IV B 1. In Fig. 5 and in Fig. 7 these 
phases arc shown as function of the regeneration rate a. 
In Fig. 8 the initial concentration of tubulin cq was cho- 
sen very close to the threshold cqc with e = 0.01. At this 
value of the control parameter the oscillations behave 
harmonically and the agreement between the numerical 
solution and the amplitude approximation is rather good, 
as described in the following section. For larger values 
of the reduced control parameter the oscillations become 
anharmonic. 
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FIG. 9. For model II the time-dependence of the polymer- 
ized tubulin r/XL (solid), ct (dashed) and Coh (dotted) arc 
shown below the Hopf bifurcation (sec Fig. 6) for a = 0.0036 
in a) and for a = 0.08 in b) with an initial concentration 
Co = 70. The corresponding critical concentration for both 
values of a is coc = 134.15. The stationary value of the tubu- 



lin-d concentration is c^°^ = 2.465 in a) and cjj"' = 34.305 in 
b). In c) the time dependence is shown beyond the oscillation 
threshold at a = 0.01. In part d) the length distribution of 
the growing microtubules Pg{l) is plotted at three different 
times. Both in c) and d) the initial concentration is Co = 80 
and the critical tubuhn concentration is coc = 65.97 which 
corresponds to the value e — 0.21 for the reduced control pa- 
rameter. The maximal length of the growing microtubules 
has been chosen as 12 • Vgj ffj^. In all parts the rest of the 
parameters axe v = 0.01, x = 0.01, Vg — 0.1. 



As already been mentioned in the introduction, the 
length distribution of filaments is a crucial difference be- 
tween the biochemical reaction discussed in this work and 



14 



the common oscillatory chemical reactions. For model I 
we show in Fig. 8d) at two different times and at e = 0.01 
the superposition of the length distribution of growing 
and shrinking microtubules, cf. P{1) = Pg{l) + Ps{l)- 
The exponential decay of the envelope of the length dis- 
tribution is described by Eq. (37b) and Eq. (37c), with 
a decay rate Vg/fj,^l, and the modulation is due to the 
traveling waves in the time-dependent contribution. The 
amplitude of the shrinking microtubules is rather small 
for /3 = 0.1, cf. Eq. (37c), and therefore the contribution 
to P{1) comes mainly from growing microtubules. 

For model II a discretization of the length coordinate 
in the equation for growing microtubule, cf. Eq. (2), 
also provides a stable numerical algorithm. Hence, the 
nonlinear oscillatory solution can be obtained numeri- 
cally without the approximations as described for model 
I in Eqs. (50) above. The respective results are shown in 
Fig. 9a)-c), where the densities i, q and Cqu are shown as 
function of time for three different regeneration rates a. 
For both values of a in part a) and b) the tubulin con- 
centration is smaller than the corresponding threshold 
value but the absolute distance cq — cqc to the thresh- 
old is equal. These transient subthreshold-oscillations 
are remarkable, because the transient oscillations in ex- 
periments might be subthreshold ones in contrast to the 
common interpretation that the oscillations are transient 
due to the tubulin-t consumption during the microtubule 
polymerization. 

For the simulations shown in Fig. 9 a narrow length 
distribution pg has been used as initial condition. In 
such cases the tubulin-t concentration corresponds al- 
most to the total initial concentration cq. Starting with 
such an initial condition, at first tubulin-t dimers are 
consumed during the growth of microtubules. This leads 
to a first maximum of the polymerized tubulin L, but 
the oligomers, the decay product of the microtubules, are 
negligible and as consequence the densities of tubulin-d 
and tubulin-t drop down too. But a small increases 
the catastrophe rate feat and microtubules decay with a 
higher rate, which increases the density of oligomers etc.. 
After a few of such oscillations the densities reach their 
stationary values for subthreshold concentrations. In the 
case of a large regeneration rate the oscillation frequency 
is much higher than for small regeneration rates and more 
oscillations are performed until the stationary values are 
reached. The stationary value of the polymerized tubulin 
and of the oligomers is also much larger in part b) than 
in part a) whereas the stationary value for the tubulin-t 
dimers remains nearly constant. The reason is the low 
density Cd in the case of large values of a. 

Far beyond the threshold of the Hopf bifurcation, the 
oscillations become anharmonic as shown in Fig. 9c). 
Hereby the oscillations of the total polymerization and of 
their decay product, CoH, are nearly in antiphase, similar 
as in experiments described in Ref. [7]. At the threshold 
a phase difference of tt was only possible in the limit of 
large regeneration rates a and for a much smaller disso- 



ciation rate x- If we consider initial concentrations which 
are much larger than the corresponding critical concen- 
tration, the phase shift of tt between L(t) and Cou{t) is 
also possible at intermediate values of x a-nd a. 

In Fig. 9d) the distribution for growing microtubules 
Pg{l) is shown for model II at three different times. The 
reduced control parameter e = 0.21 is rather large com- 
pared to its value in Fig. 8 and the curves indicate that 
the length distribution of microtubules isn't described 
anymore by harmonic traveling waves as in the vicinity 
of the bifurcation. As long as the tubulin-t density is 
large and the catastrophe rate small, microtubules grow 
with a constant velocity Vg and only a few of them ex- 
perience a catastrophe. During this period a plateau in 
the length distribution is build. But after a large amount 
of tubulin-t has been used up the catastrophe rate feat 
increases very steeply leading to a strong decay of the 
microtubules at all lengths. If fcai drops down again the 
low density of long microtubules grow further with a con- 
stant velocity Vg and short microtubules are nucleated at 
a higher density. This leads to a step in the distribution 
that travels with the growth velocity Vg to larger values 
of I. Therefore, far beyond the oscillatory threshold the 
temporally anharmonic behavior of Ct leads in this man- 
ner to step-like length distribution of the microtubules. 



VI. AMPLITUDE EXPANSION 

Here we focus on a semi-analytical treatment of the 
oscillating polymerization slightly beyond its onset. The 
interesting question here is whether the bifurcation to 
these oscillations is continuous (supercritical) or discon- 
tinuous (subcritical). In physical systems Hopf bifurca- 
tions are mostly subcritical [27,35], but for the models 
discussed in this work we always find a supercritical one. 
This is advantageous, because for a supercritical Hopf 
bifurcation a semi-analytical treatment is possible. The 
appropriate frame work is the universal amplitude equa- 
tion of oscillatory fields, cf. Eq. (1). This equation will 
be derived in the present section from the basic reaction 
equations of microtubule polymerization. 

The perturbation analysis employed for the derivation 
of Eq. (1) is an expansion of the solutions of the basic 
equations with respect to small amplitudes of the oscil- 
latory contributions [27,28]. As a small perturbation pa- 
rameter the relative difference between the actual tubulin 
concentration cq and the critical tubulin concentration 
coc is introduced 



COc 

A signature for the ± symmetry of the oscillatory behav- 
ior [27,28] is the power law for the oscillation amplitude 

A ^ ^fe. Accordingly the solutions of the basic equations 
at the threshold are expanded with respect to powers of 
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U = + eV2u(l) + ^ u(2) + ^3/2^(3) ^ ^(^2^ ^ (gg) 



roat^ = s{l + ia)A- g{l + ic) \A\'^A 



(64) 



di) 



is used 



where the vector notation u^-'' = {pg ,p. 
with j = 0, 1, 2, 3. The components of u^^^ differ by a fac- 



tor of ^/e, such as ^/s c[^^ = cY' etc. . The components of 
u^°-* describe the stationary microtubule polymerization 
as given in Sec. Ill and the components of u^^^ describe 
the linear oscillatory solutions that may be written at the 
threshold in the following form 



(1) 



(61) 



Here Ue includes the amplitude ratios between the fields 
c\^\ p^g \ p^P at the threshold and \feB = A as ex- 
plained below. 

Close to the threshold one has i?e(a) ~ e <C 1 and 
the linear solution u^^^ ^ e°'* grows or decays only by a 
very small amount during one oscillation period 2'k/ijJc- 
These two disparate time scales near the threshold, that 
of the oscillation period (oc 2tt/lOc) and that of growth 
and decay (oc 1/e), may be separated within a perturba- 
tion expansion by introducing a slow time scale T = et 
[27,28]. The fast time scale is included in the exponential 
function e*'^'=* and the other one will be described by a 
slowly varying amplitude B{T). Accordingly the linear 
solution near threshold may be written as 



u(i)(t,T) = B(r) Ue e*'"=*-hc.c. 



(62) 



In order to differentiate this product of time dependent 
functions, instead of applying the chain rule of differen- 
tiation, one may replace this operation by the following 
sum dt dt + edr- Here dt acts only on the fast time- 
dependence occurring in the exponential function and Ot 
acts only on the amplitude B{T). 

Using this replacement and the e-expansion of u the 
basic equations given in Sec. II can be ordered with re- 
spect to powers of \/e. In this way we obtain a hierar- 
chy of partial differential equations, which we need up to 
0{£^/'^). The whole procedure is described in more detail 
in Appendix A. The amplitude equation follows from a 
solubility condition for the equation at order 0{£^^'^ ) and 
it has the following form 



TodrB ^ {l + ia)B - g{\ + ic) \B\'^B . 



(63) 



To is the relaxation time, a is the linear and c is the non- 
linear frequency shift. The nonlinear coefficient g deter- 
mines the bifurcation structure. For .9 > the bifurcation 
is supercritical (steady) and for g < the bifurcation is 
subcritical (unsteady). For the coefficients tq, a, g and c 
one obtains long expressions in terms of the reaction con- 
stants of the basic equations, that have been calculated 
by using computer algebra. The respective formulas are 
not presented but the parameter dependence of the coef- 
ficients is shown in Fig. 10 for model I and in Fig. 11 for 
model II. 

Rescaling the time T = et and amplitude A = ^/e B 
back to the original units yields the amplitude equation 



as introduced in Sec. I. This equation has simple nonlin- 
ear oscillatory solutions of the form A = f e*^* , with an 
amplitude F and a frequency f2 as follows 



F = 



n ^ — {ea- gcF'^) 



TO 



(65) 



il describes the deviation of the oscillation frequency 
from the critical one, u)c- 

The linear coefficients tq and a of Eq. (64) may directly 
be calculated from the dispersion relation in Eq. (29) or 
Eq. (44) in the following way. The solution ^ = of 
Eq. (64) corresponds to the stationary polymerization 
described in Sec. HI, which is in the range £ < stable 
against small perturbations A ~ Fe"^* ( with F <^ \/|e|) 
and unstable for e > 0. Neglecting in Eq. (64) the con- 
tributions due to the cubic nonlinearity one obtains from 
its linear part the dispersion relation a = e{l + zo)/ro. 
This formula gives the relaxation time tq and the lin- 
ear frequency dispersion a in terms of derivatives of the 
growth rate with respect to the control parameter e: 
To = {dRe{(7)/de)~^ and a = TodRe{a) / de. If e is ex- 
pressed in terms of the dimer density co, cf. Eq. (59), 
then both quantities may also be written in terms of the 
derivatives with respect to cq 



1 



cocdRe{a)/dco 



COcTO 



dlm{a) 
dco 



(66) 



whereby both derivatives are taken at the threshold con- 
centration cqc- The dispersion relation o"(e) as obtained 
on the one hand by the amplitude equation and on the 
other hand by solving Eq. (29) or Eq. (44), both have 
to reproduce the growth or decay dynamics of small per- 
turbations with respect to the stationary polymerization. 
Therefore, the coefficients tq and a of the amplitude equa- 
tion can directly be calculated via Eq. (66) from the nu- 
merical solutions of Eq. (29) or Eq. (44). 

One aim of the amplitude expansion is the determina- 
tion of the type of the Hopf bifurcation. For two different 
nucleation rates i' = 0.01 and = 0.04 the variation of 
the nonlinear coefficient g as function of the regenera- 
tion rate a is shown for model I in Fig. 10 (top) and 
for model II with oligomer dynamics in Fig. 11. In both 
cases g behaves rather similar and g is positive for the 
models investigated in this work. Therefore the Hopf 
bifurcation is supercritical. In Fig. 10 and Fig. 11 the 
nonlinear coefficient g increases at first with the regener- 
ation rate a and reaches a maximum in a range where the 
threshold concentration coc{a) takes its minimum. The 
corresponding threshold curves cqc{oi) for two different 
nucleation rates ly = 0.01 and ly = 0.04 also cross each 
other. Beyond this maximum of g { where coc{a) takes 
its minimum) decreases again. 
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It should be mentioned that for a given vahie of the 
control parameter e a large value of g corresponds to a 
small value of the oscillation amplitude. Since the thresh- 
old coc(q;) varies too, the variation of the oscillation am- 
plitude with a is much less when cqcE is kept constant. 
According to the sign of the nonlinear coefRcient c the os- 
cillation frequency Wc+f^ decreases with increasing values 
of e. The results in terms of the amplitude equation are 
in fairly good agreement with the behavior of the full 
numerical solution of the basic reaction equations. 
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FIG. 10. The coefficients to, a, g and c of the amplitude 
equation (64) are shown for model I as function of the regen- 
eration rate a and for two different nucleation rates v = 0.01 

(solid line) and p — 0.04 (dashed). For the rest of parameters 
the values Vg = 0.1, /? = 0.1, / = 0.1 and c/ = 3 have been 
chosen. 



A determination of the bifurcation structure by nu- 
merical simulations of the basic equations is error prone 
compared to results of perturbation calculation described 
here. Besides the lower accuracy, parameter studies such 
as in Fig. 10 and Fig. 11 are with numerical simulations 
much more time consumptive. 
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FIG. 11. The linear and nonlinear coefficients of the am- 
plitude equation (64) are shown for model II as a function 

of the regeneration rate a and for two different nucleation 
rates u = 0.01 (solid) and v = 0.04 (dashed). For the rest of 
parameters the values x = 0.01, / = 0.1, c/ = 3 have been 
chosen. 

Close to threshold the advantages of the perturbation 
calculation are obvious. However, it is a priori unknown 
in which e-range the amplitude equation approach (64) 
applies quantitatively. For some systems the amplitude 
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equation is valid in a rather large range of the control pa- 
rameter £, but for other systems its validity is restricted 
to very small values of it, cf. Ref. [27]. In order to check 
this for our models of microtubule polymerization, we 
compare in Fig. 12 the variation of the oscillation am- 
plitude of Cj^^ with the control parameter e as obtained 
by the numerical solution described in Sec. V and by the 
solution A = sj g of the amplitude equation for two 
different values of the nucleation rate v. At larger values 
of the control parameter, e = 0.1, the difference between 
the results for the numerical solution and the amplitude 
equation is still less than 8%. For model II the devi- 
ations are larger between the amplitude determined by 
the amplitude equations approach and the numerical so- 
lutions with ansatz in Eq. (52). However, when we solve 
the equation for growing microtubule, cf. Eq. (2), nu- 
merically by discretization of the length coordinate, the 
deviations becomes smaller. 




control parameter e 

FIG. 12. The amplitude of the oscillations as function of 
the reduced control parameter e and for two different nucle- 
ation rates v. The solid line is the result of the amplitude 
equation and data points have been determined from simu- 
lations as described in section V. The parameters that have 
been used are 7 = 1, a = 0.01, Vg = 0.1, f3 = 0.1, / = 0.1, 
Cf = 3. 

For both models, the linear coefficient ro and the non- 
linear coefficients g, c don't differ very much from each 
other. However, the linear frequency shift a increases in 
model II for large regeneration rates a whereas for model 
I it decreases. Nevertheless in both models the nonlinear 
frequency correction due to the values of c is much larger 
than the linear correction due to a, c >> a. 

Whether the bifurcation to oscillatory polymerization 
is also supercritical in experiments under stationary re- 
generation conditions, is an open question. Therefore 
an experimental determination of the bifurcation-type 
would be an important test for the reduced models in- 
vestigated in this work. 

The variation of the other coefficients with the re- 
generate a provides further contact between the model 



parameters and experimentally measurable quantities. 
The parameters c, a and tq may be determined as 
follows. Studying the growth of small perturbations, 

cj.^^ cx e-^-'^^'^^* = e"^*/"^" by plotting the logarithm ef/ro cx 

log(cj^') as function of time and for different values of e, 
the relaxation time to may be determined. In a similar 
manner a may be determined by studying the frequency 
of a perturbation far below its nonlinear saturation am- 
plitude. If the perturbation saturates finally, the oscilla- 
tion frequency of the nonlinear solution changes with e as 
indicated by Eq. (65). From this e-dependence the non- 
linear coefficient c may be extracted. An experimental 
determination of these coefficients, as described, would 
be a further test of the basic model equations. 



VII. SUMMARY AND CONCLUSION 

Two reduced models, that capture oscillating micro- 
tubule polymerization and the length distribution of the 
microtubule filaments, as described in Sec. II, have been 
analyzed. In both models the complex biochemical re- 
action steps of microtubule polymerization are described 
by a few ones, which have been identified in experiments 
to be important. The focus on a few essential degrees of 
freedom leads to some simplicity of the models that al- 
lows for instance a derivation of analytical expressions for 
the threshold and the oscillation frequency at the Hopf 
bifurcation. Such analytical results make trends as func- 
tion of the reaction rates more easily visible. Some of 
these trends may be tested in experiments and some of 
the reaction constants may be measured. 

At threshold also analytical expressions could be de- 
rived for the temporal evolution of the concentrations 
and the length distribution of microtubules. Those pro- 
vide a detailed picture about the temporal variation of 
the fields, their relative phases and the amplitude ratios 
between them. The formula for the length distribution 
is especially instructive, cf. Eq. (37b), it describes a su- 
perposition of amplitude oscillations of the distribution 
and traveling waves, where the waves always travel to- 
wards larger lengths. This qualitative behavior of the 
length distribution during oscillatory polymerization is 
rather independent of the respective model and it is a 
rather general feature. A few snap shots of the numer- 
ically generated distribution far beyond threshold are 
shown in Fig. 9b). The distribution includes still travel- 
ing waves, but far beyond threshold these behave rather 
anharmonic. 

For stationary reaction conditions, as assumed in this 
work. Fig. 9 shows a remarkable result. In part a) and 
b) of this figure a subthreshold concentration for tubu- 
lin was assumed, i.e. cq < cqc, and in both cases the 
final state is stationary polymerization. However, on the 
route to the stationary state transient oscillations occur. 
Therefore, the transient character of the microtubule os- 
cillations observed in an experiment with an enzymatic 
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regeneration process for GTP [11], could have, according 
to the results described in this work, its origin in a to low 
tubulin concentration. A higher tubulin concentration or 
an appropriate regeneration rate of GTP and a different 
live time of oligomers could lead in a similar experiment 
to persistent microtubule oscillations. 

For transient oscillations observed in other experi- 
ments the common interpretation is as follows. During 
the microtubule polymerization in experiments the avail- 
able GTP is used up and the oscillations last only for 
a few periods. If GTP is continuously supplied, the si- 
multaneously increasing amount of GTD inhibits various 
reactions steps and slows down the reaction cycle. The 
results shown in Fig. 9 a) and Fig. 9 b) indicate that os- 
cillations may occur as a transient because either GTP 
is used up or the tubulin concentration has only a sub- 
threshold value. Accordingly, there may be several rea- 
sons for transient oscillations in experiments. Either the 
initial tubulin concentration has a subthreshold value, 
the decay rate of oligomers and the regeneration rate for 
GTP do not have their optimal values which would ex- 
plain transient oscillations in experiments with a regen- 
erative enzyme system, the reaction conditions arc not 
constant, because GTP is used up. In the latter case the 
available tubulin-t decreases with time and the micro- 
tubule polymerization decays. 

In an in vitro experiment with constant reaction con- 
ditions the threshold of oscillations might be measured 
by increasing the tubulin dimer concentration by appro- 
priate steps. Immediately after each step transient os- 
cillations might occur, but the threshold is only crossed 
when the oscillations persist over a long time. 

In order to avoid numerical instabilities during long 
time simulations of the reaction equations, including 
Eq. (6b) , we use analytical approximations for the length 
dependence of the microtubule distributions as described 
in Sec. V. This stable numerical scheme can be general- 
ized in future work to an effective algorithm for dealing 
with microtubule polymerization in one and two spatial 
dimensions [26] in order to investigate spatial patterns 
occurring during polymerization of microtubules [10,24]. 
The respective extension of the amplitude equation may 
also lead to new interesting insights. 

Microtubule filaments at a high density show a 
isotropic-nematic phase transition [22], similar as it has 
been observed for F-actin filaments [36]. Since the 
early theory of Onsager [37] this transition for rods in 
a solvent is a well understood phenomenon [38]. For a 
monodisperse distribution of filaments of fixed length, 
i.e. without polymerization kinetics, many aspects of the 
isotropic-nematic transition have been understood. For 
polydisperse rod mixtures some aspects of the isotropic- 
nematic transition can be considered to be understood 
too [39]. The effects of nucleation, growth of filaments 
and the decay of filaments on the isotropic-nematic 
transition are not known at present and one may ex- 
pect interesting phenomena related to this kinetics [23]. 
The effect of oscillating microtubule polymerization on 



the isotropic- nematic transition is also completely unex- 
plored at present and will be investigated in forthcoming 
works [23]. 

It is a great pleasure to thank M. Breidenich, H. Fly- 
vbjerg, E. Mandelkow, H. MiiUer-Krumbhaar, J. Prost 
and J. Tabony for interesting discussions. 



APPENDIX A: AMPLITUDE EXPANSION FOR 
MODEL I 



For model I and the catastrophe rate given in Eq. (3) 
the major steps of the derivation of the amplitude equa- 
tion (1) arc described in this appendix. Since we neglect 
the rescue of shrinking microtubules, fresc = 0, the only 
nonlinear term in the basic equations of model I is the 
product fcat{ct)Pg in Eq. (6). At first we expand the 
concentrations Ct,co and length distributions Pg^s with 
respect to deviations from their stationary value at the 
threshold for oscillation, such as for instance for the 

tubulin-t concentration, i.e. Ct — c|°^ = ^c[^^ +ec[^^ + 

Note that the tilded fields differ just by a power of 
from the untilded fields as introduced in Sec. IV, cf. 



„(i) 



etc. . In order to simplify the notation 



of this appendix we drop the "tilde" and the expansion 
of the catastrophe rate takes the form 

fca. = fn+ e'"f^^ + efT^ + e^"f^. + ■•■, (M) 
whereby the coefficients of this expansion are 



/■(I) 

J cat 



/■(2) 

J cat 



f(3) 
J cat 




AO) 

J cat 



(A2a) 
(A2b) 



— - ^ — • (A2c) 



Collecting in Eq. (6) and Eq. (11) the contributions to 
the order e^^"^ we recover the linear equations given in 
Sec. IV 

dtPi'^ = f'^t - /ia] P« - vM^^ , (A3a) 

dtPi'^ = + /iS + vsdw^^^ , (A3b) 



dv 



{' / \ 

cW=-7 / ci^(^g?>W+a«(p«+pW)) 



(1) 



(A3c) 



At order e we obtain the three equations 
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' (2) 

f(0)f^„(0) _ f(").„(2) 
J cat -Pg ./<;a( -f-^g 



'J cat „ Pg 
Cf 

/ (2) 

-/■(0)£^„(0) 
J cat „ -r^g 



(1) J-(O) , 




+acoc 



ac. 



(2) 



(A4a) 



(A4b) 



(A4c) 



—ac. 



(3) 



(A6c) 



The two fields p^^ and p^^^ have to be calculated ex- 
plicitly at this order from Eq. (A6a) and Eq. (A6b) 
as well. With both solutions the integral on the right 
hand side of Eq. (A6c) can be calculated. Eq. (A6a) 
and Eq. (A6b) include both contributions proportional to 



and 



but only the single harmonic terms are 



relevant in Eq. (A6c). The coefficient of e*'^'=* in Eq. (A6c) 
must vanish. Part of it vanishes automatically, because it 
reproduces the threshold condition and the rest provides 
the amplitude equation with all the coefficients now given 
in terms of the reaction rates of the basic equations. 



With the solutions of the equations at the previous order 

e-*^/^, that arc already given in Sec. IV A 3, the equations 
at order e have to be solved. These solutions have the 
following form 



„(2) 



Aq + A2 cxp {2iuJct) + c.c. 



P. 



.2iuJct 



(A5a) 

g - ^ v""vv ' v-^y-r- ■ cc] ) , (A5b) 

(2) ^ ^-fi:il/v, (^jy^^i^ ^ [i?2(/)e2*'^^* + C.C.]) , (A5c) 



whereby the expressions for the coefficients Aq, A2, 
and Di are rather lengthy in terms of the coefhcients of 
the solutions at order e^^"^ and are not given here. 
The equations at the next higher order e"^/^ are 



(3) 
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